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Abstract. The structure of steady plane-parallel radia- 
tive shock waves propagating through the hydrogen gas 
undergoing partial ionization and excitation of bound ato- 
mic states is investigated in terms of the self-consistent so- 
lution of the equations of fluid dynamics, radiation trans- 
fer and atomic kinetics. The shock wave model is repre- 
sented by a flat finite slab with no incoming radiation 
from external sources at both its boundaries. The self- 
consistent solution is obtained using the global iteration 
procedure each step of which involves (1) integration of 
the fluid dynamics and rate equations for the preshock 
and postshock regions, consecutively, both solutions being 
fitted by the Rankine-Hugoniot relations at the discontin- 
uous jump; (2) solution of the radiation transfer equa- 
tion for the whole slab. The global iteration procedure 
is shown to converge to the stable solution which allows 
for the strong coupling of the gas flow and the radiation 
field produced by this flow. Application of the method 
is demonstrated for the shock waves with upstream ve- 
locities of 15 km s _1 < Ui < 60 km s _1 (i.e. with up- 
stream Mach numbers 2.3 < M% < 9.3) and the hydrogen 
gas of unperturbed temperature T = 3000K and density 
p = 10~ 10 gm cm~ 3 . 

Key words: Shock waves - Hydrodynamics - Radiative 
Transfer - Stars: atmospheres 



1. Introduction 

Radiative shock waves belong to the conspicuous phe- 
nomena demonstrating the tight interplay between hy- 
drodynamic motions and the radiation field. The role of 
this interplay is strongest in the low gas density flows, so 
the shocks are of tremendous importance in astrophysics. 
They are observed in a wide variety of astrophysical phe- 
nomena: nova and supernova explosions, bright filaments 
in old supernova remnants, accretion flows in protostel- 
lar clouds. Shock waves are detected also in atmospheres 



W Vir, RV Tau and Mira type stars. Periodic shocks prop- 
agating through pulsating atmospheres lead to the disten- 
tion of outer atmospheric layers and to the mass loss. 

Importance of radiative shock waves attracted atten- 
tion of many authors but nevertheless the shock prop- 
erties are explained quite well still qualitatively (see, for 
example, Zel'dovich & Raizer 1966; Skalafuris 1968; Mi- 
halas & Mihalas 1984; Liberman & Velikovich 1986). The 
principal difficulty in obtaining the correct quantitative 
description of the shock wave structure is that the model 
has to allow for the strong coupling between the gas flow 
and the radiation field, both them being characterized by 
substantial departures from LTE. Solution of this prob- 
lem encounters serious difficulties, so that in immensely 
numerous studies available at present in the literature the 
authors used various assumptions and simplifications (e.g. 
local thermodynamic equilibrium, treatment of the radi- 
ation transfer in diffusion approximation, neglecting the 
opacity in the Balmer continuum etc.). Many of these as- 
sumptions were found later inadequate or leading to un- 
certain conclusions. For instance, Kogure (1962), Sachdev 
(1968) and Hill (1972) used the LTE approximation which 
does not hold as emphasized later Narita (1973). Consid- 
ering the hydrogen gas, Whitney & Skalafuris (1963) re- 
laxed this assumption but incorrectly assumed that the 
postshock region is transparent for all hydrogen continua. 
Finally, Narita (1973) took into account the opacity in 
the both Lyman and Balmer continua. The most elabo- 
rate numerical modelling of radiative shock waves based 
on the self-consistent solution of the equations of fluid 
dynamics, radiative transfer and atomic level populations 
was done by Klein et al. (1976, 1978). However, the coarse 
zonning (~ 10 6 cm) did not allow to authors to consider 
the detailed structure of the shock front including the ra- 
diative precursor and the thermalization zone where the 
electron temperature gradually equalizes with tempera- 
ture of heavy particles. Nevertheless, this approach was 
found to be enough for consideration of shock dynamics 
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cursor is not so important due to the high temperature of 
the unperturbed gas. 

Because astrophysical shocks in stellar atmospheres 
propagate through the partially ionized hydrogen g as, a 
substantial fraction of photons produced within the wake 
are absorbed in the radiative precursor. As was shown by 
Gillet & Lafon (1984, 1990) the structure of the radiative 
precursor is complex and should be treated with same de- 
gree of approximation as the postshock region. In their 
studies Gillet & Lafon (1984, 1990) treated the radiative 
transfer as an initial value problem which was solved us- 
ing the shooting method. The principal obstacle in such 
an approach is that the transfer equation possesses a sin- 
gularity in the postshock region (Gillet et al. 1989). In- 
deed, application of the eigenvalue methods for solution 
of the nongrey transfer problem is affected by exponen- 
tially growing errors (Mihalas 1978). 

One of the first attempts to obtain the self-consistent 
solution for the shock wave structure was undertaken by 
Nelson & Goulard (1969) and Nelson (1973). They con- 
sidered the shock waves propagating through the argon- 
like gas with upstream Mach numbers of M\ = 18 and 
Mi = 24. The continuity, momentum and energy equa- 
tions were written in the integral representation whereas 
the radiation transfer was treated in the simplified formu- 
lation. The studies of radiative shock waves in helium and 
nitrogen done by Clarke & Ferrari (1965), Farnsworth & 
Clarke (1971) and Foley & Clarke (1973) seem to be the 
best among known in the literature. The authors empha- 
sized the crucial role of the radiation transfer treatment 
and employed the formal solution of the transfer equa- 
tion. The self-consistent shock wave models were obtained 
in these studies with iteration procedure. Unfortunately, 
there is a problem of exponential factors when the formal 
solution is applied for optically thick layers. 

In this paper we present a new approach based on 
the iterative solution of the equations of fluid dynamics, 
the rate equations and the radiation transfer equation. 
The momentum equation, the energy and rate equations 
are written in the form of ordinary differential equations. 
These equations are stiff and such a representation is most 
appropriate from the point of view of stability and small 
truncation errors. The radiation transfer is treated as a 
two-point boundary value problem. This allows us to ob- 
tain the stable solution of the transfer equation for the 
whole spectral range including both the opaque Lyman 
continuum and the more transparent higher order con- 
tinua. The method of global iterations takes into account 
the coupling between the gas flow and the radiation field, 
so that the structure of the radiative shock wave is con- 
sidered in terms of the self-consistent model. 

In the framework of this first approach, only devoted to 
provide a new technique for obtaining the self-consistent 
solution, we consider the structure of steady, plane-parallel 
shock waves propagating through an infinite, isotropic, 



for example, to a good degree of approximation in most 
applications to stellar atmospheres. Indeed, the time re- 
quired for the gas flow to cross the shock wake with typical 
thickness of 10 4 -10 7 cm is much less than the characteristic 
time during of which the shock wave energy appreciably 
decreases. In pulsating stars, the radiative lifetime of a 
shock is between a few hours and a few months, which, 
consequently, is much larger than the 10~ 2 -10 s of the gas 
flow time to cross the shock wake. Accuracy of the plane- 
parallel approximation follows from the very small width 
of the shock wave in comparison with stellar radius. 

After describing the shock wave model (Sect. 2) we de- 
rive the system of ordinary differential equations (Sect. 3). 
The radiation transfer equation is solved for the whole 
shock wave model using the Feautrier technique (Sect. 4). 
In Sect. 5 we show that the global iteration procedure com- 
prising the initial value problem for ordinary differential 
equations and the two-point boundary value problem for 
radiation transfer converges to the self-consistent solution. 
Results of calculations demonstrating the applicability of 
the method are given in Sect. 6. Finally, in Sect. 7, we give 
some concluding remarks and discuss the future aspects 
of the problem. 

2. The shock wave model 

Consider a steady, plane-parallel shock wave propagating 
through the homogeneous medium which is at rest and 
consists of a pure hydrogen gas. No radiation and grav- 
itational forces from external sources are assumed to be 
present. The problem to be solved is that to describe the 
spatial structure of the shock wave in terms of the self- 
consistent solution of the equations of fluid dynamics, the 
rate equations for hydrogen atomic level populations and 
the radiation transfer equation. The problem is character- 
ized by three input parameters: the temperature T\ and 
the density pi of the unperturbed gas as well as the speed 
Ui at which the gas material flows into the shock. 

The structure of radiative shock waves is schemati- 
cally divided into four zones: (1) a precursor, where the 
gas is heated and is partially ionized by radiation emerg- 
ing from the postshock region; (2) a very narrow zone, 
where a major part of the kinetic energy of the upstream 
flow is converted due to viscosity and conductivity into 
the thermal energy of translational motions of heavy par- 
ticles, that is, neutral atoms and ions; (3) a thermalization 
zone, where the kinetic energy of translational motions of 
heavy particles is redistributed among various degrees of 
freedom; (4) a radiative relaxation zone, where hydrogen 
atoms rccombine and the gas radiatively cools. Because of 
its extremely small width (a few mean free paths of gas 
particles), the second zone cannot be correctly described 
in terms of the fluid dynamics and, hence, should be con- 
sidered as a discontinuous jump across which the Rankine- 
Hugoniot equations are applied. Thus, the present study 
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radiative precursor as well as the thermalization and re- 
laxation zones. 

Let the origin of the comoving frame to coincide with 
infinitesimaly thin discontinuous jump dividing the me- 
dium into the preshock and the postshock regions. The 
spatial coordinate is X — at the discontinuous jump, 
is negative in the preshock region and is positive in the 
postshock region. The shock wave model is represented by 
a flat finite slab comoving with discontinuous jump. Thus, 
the velocity of the gas material flowing through the slab 
is always positive: U = dX/dt > 0. At the outer boundary 
of the preshock region with spatial coordinate X\ < the 
gas is assumed to be unperturbed. The spatial coordinate 
of the postshock outer boundary is X^ > and physical 
properties of the gas at this point are not known. 

The radiation transfer equation is solved for the whole 
slab in the framework of the two-point boundary value 
problem, therefore the slab is represented by a set of spa- 
tial cells. The discontinuous jump locates at the J-th cell 
boundary with spatial coordinate Xj = 0. The cells are 
incqually spaced, are smallest at the discontinuous jump 
and increase outwardly in both directions from the dis- 
continuous jump according to the geometrical progression. 
Following the traditional conventions of computational ra- 
diation hydrodynamics (see, for example, Mihalas & Mi- 
halas 1984), all thermodynamic variables are defined at 
the cell centers Xj_i/ 2 — \ {Xj-\ + Xj) and are denoted 
by half-integer subscripts. The spatial coordinates of the 
cell centers nearest to the discontinuous jump hereafter 
are referred as the inner boundaries of the preshock and 
postshock regions, respectively, and are denoted as Xj_ 1 / 2 
and X J+1/2 . 

The assumption of the steady shock wave allows us 
to reduce the equations of fluid dynamics to a system of 
ordinary differential equations. In order to take into ac- 
count the coupling between hydrodynamic motions and 
the radiation field we employ the iteration procedure com- 
prising the consecutive solution of the transfer equation 
and integration of the rate and fluid dynamics equations. 
The starting point for integration of ordinary differential 
equations is the outer boundary of the preshock region X\ 
where the gas is assumed to be unperturbed. The preshock 
integration is done within the interval [A"i, ATj_ 1 / 2 ]- Then 
we solve the Rankine-Hugoniot equations 

gll = m = Co , 
mU + P g = Ci 
'h 1 
.Q 2 
where h is the enthalpy, 

h = 5/3E tr ans + E rot + E v u, + E e \ + Ediss + Ei, 

where E trans , E roU E vib , E eh E diss and E io „ are the trans- 
lational energy; rotational, vibrational, electronic excita- 
tion internal energies and dissociation and ionization po- 
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(2) 
(3) 
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Ei on are different of zero. F T is the total radiative flux, Co, 
Ci and C 2 are the mass, momentum and energy fluxes 
across the discontinuous jump. In the present study we 
assume that the radiation pressure P ra d and the radia- 
tion energy density £7 ra d can be neglected in comparison 
with gas pressure and internal energy of the gas material. 
Applicability of this assumption is shown below. 

Let us designate, for the sake of convenience, the quan- 
tities defined at A^j_!/ 2 by the superscript minus and 
the quantities defined at Xj + i/ 2 by the superscript plus. 
Eqs. (|l|) - (g) are solved for the temperature of heavy 
particles T+ and the inverse compression ratio lo = p~ / p + 
= U + /U~. Assuming that hydrogen atoms undergo across 
the discontinuous jump neither excitation of the bound 
levels nor ionization and that the hydrogen ions have the 
same temperature as that of neutral hydrogen atoms, we 
obtain 



T~ 



2F; 



t: 



'H 

-F+ 



T c + ) 



1 rhU , 
5 n^fc 



5 U n H k 



(5) 



where k is the Boltzmann constant, ?t,h and n e are the 
total numbers of hydrogen atoms and free electrons per 
unit volume, respectively. It should be noted that jump 
conditions (1) - (3) are not applied to the infinitasimally 
thin discontinuous jump coinciding with the cell interface 
X] but rather relate the physical variables at cell cen- 
ters Xj_i/2 and Xj +1 / 2 . That is why the term containing 
F r ~ — F T + into Eq. (5) does not cancel and has to be taken 
into account. The inverse compression ratio is obtained 
from 



Au? + Bw + C = , 



where 



.4 



-mil 



B = -rhZT - P t 
1 
5 

Eqs 



C 



-mil 



2F± 
5 



(6) 

(7) 
(8) 
(9) 



and (J6|) imply that the postshock electron tem- 
perature T e + at the cell center Xj +1 / 2 is known. Across the 
discontinuous jump the elctron gas undergoes the adia- 
batic compression and the electron temperature increases 
by a factor of where 7 is the ratio of specific 

heats (Zel'dovich & Raizer 1966). Because the adiabatic 
compression has a weak effect and the electron heat con- 
duction is out the scope of the present study, we assumed 
that the electron temperature does not change across the 
discontinuous jump, that is, T~ = T c + . 

Eqs. (g) and (g) are solved to determine the initial con- 
ditions for the postshock integration which is done within 
the interval [Xj + i/ 2 , X^]. Integration of ordinary differ- 
ential equations provides with spatial distibutions of elec- 
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free electrons n c and atomic level populations m. These 
quantities are used in solution of the radiation transfer 
equation which gives improved radiation intensities within 
the whole slab. The global iteration procedure consisting 
of the integration of ordinary differential equations and 
solution of the radiation transfer equation continues un- 
til the relative changes of all quantities become less than 
the convergence criterion. Below we discuss details of the 
global iteration procedure. 

3. The system of ordinary differential equations 

The solution vector to be found from integration of the 
system of ordinary differential equations consists of L + 4 
variables and is 



Y = {U, E' a , E' e , n' c , n[. 



(10) 



where L is the number of bound atomic states treated in 
non-LTE, U is the gas material velocity with respect to 
the discontinuous jump, 

E a _ 3 njjkTg 
~ 2 



E' 



(11) 



pip 

is the specific translational kinetic energy of heavy parti- 
cles (neutral atoms and ions), 



! E c 3 n c kT c 



(12) 



is the specific translational kinetic energy of free electrons, 
n' c is the number of free electrons per unit mass, is the 
number of hydrogen atoms in the i-th state per unit mass. 
Hereafter the prime implies that the quantity is expressed 
per unit mass. The solution vector Y does not contain the 
gas density p because this variable can be easily evaluated 
from the mass conservation relation ([]]). 

Thus, the system of ordinary differential equations 
consists of the momentum equation, two energy equations 
for heavy particles and free electrons, the rate equation 
for free electrons and L rate equations for non-LTE bound 
states of the hydrogen atom: 

dU ldP e 



dt 



p dX 



dt ~ a dt Wolc - 

dt ~~ ° dt + ^ olc ' 



dn' c 
dt 



1 L 



i=l 



dn': 1 

— = - {n e P M - riiPik) , 
dt p 



(i = l,...,L) 



(14) 
(15) 

(16) 
(17) 



where P a — nnkT a is the partial pressure of heavy par- 
ticles, P c — n c kT c is the electron pressure, Qeic and Q lnc 
are the rates of energy gain by electrons in elastic and in- 



that is collisional plus radiative, ionization and recombi- 
nation rates, respectively. It should be noted that in the 
present study we consider only bound-free transitions, so 
that equations ( |l6[) and ( |l~7|) contain only ionization and 
recombination terms. 

The system of ordinary differential equations ( |i~3| ) - 
( |l~7j ) written in the form of derivatives with respect to time 
t is not appropriate for calculation of the spatial structure, 
so that these equations should be rewritten in order their 
left-hand sides are replaced by derivatives with respect to 
the spatial coordinate X. Furthermore, the space deriva- 
tive of the gas pressure P g in Eq. (|T^ ) and the time deriva- 
tive of the specific volume V in Eqs. jlj] ) and (|l5|) have to 
be expressed in terms of integrated variables. To this end 
we write the gas pressure as a sum of translational kinetic 
energies: 



P R = nukT a + n e kT e 



-E n 



:E R 



(18) 



whereas the time derivative of the specific volume is de- 
termined from Eq. (|2|): 



dV_ 
~dt 



J_dPz 

m, 2 dt 



(19) 



Expressing the time derivative of the specific volume in 
terms of the gas pressure according to Eq. (19), P g ac- 
cording to Eq. (18) in terms of integrated quantities E' a 
and E' c and substituting Eqs. (18) and (19) into Eqs. (13) 
- (17), we obtain the following system of ordinary dif- 
ferential equations with right-hand sides depending only 
on the independent variable X and integrated variables 
{U, E' a , E' c , n' c , n[, n' L }: 



dU A 
dX ~ ~U 

^ = B (1-C ) 
dX ~ ^ aLc 



dn' e 
dX 

<H _ 
dX 

where 
A = 



B c ) , 
B C 

B c (i - c a ) , 

n e Pki 



i=l 

f^cPki T^iPik 

rh 



(z = !,...,£) 



U 2 



3U 2 



B tl 



1-P 2 



B R = 



U 1 

Qelc 



IP 2 



(20) 
(21) 
(22) 

(23) 
(24) 

(25) 
(26) 



XH 



VF r \ 



dn' c 
dX 

1 



L 

E 



i-- 



dX 



(07\ 
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C a = 
Co = 



APg 

Urn ' 
AP C 
Urn ' 



(28) 
(29) 



/3 = clt/U, cit — y/Pg/p is isothermal sound speed, xh = 
13.598 eV is the ionization potential of the hydrogen atom. 

In obtaining Eq. (|27|) we expressed the rate of energy 
gain by electrons in inelastic collisions as (Murty 1971) 



at p 



(30) 



where Ej and E' eyi are ionization and excitation energies 
per unit mass, 



V-F r = 47r J (r] u - k v J v ) dv 



(31) 



is the divergence of radiative flux, J v is the mean inten- 
sity of radiation, r\ v and k v are the total emission and 
absorption coefficients. 

Free electrons acquire the energy from heavy particles 
in elastic collisions with hydrogen ions and neutral hydro- 
gen atoms, hence, 



clc — Qci Qc 



(32) 



where Q c \ and Q oa are the corresponding rates of energy 
gain. The rate of energy gain by electrons in elastic colli- 
sions with hydrogen ions per unit mass is (Spitzer & Harm 
1953) 



2 n e T a T c 



3 Q teq 

where t cq is the time of equipartition given by 



252TL 



3/2 



n e In A 

and 

In A = 9.43 + 1.151og (T c 3 /n c ) 



(34) 



(35) 



The rate of energy gain by electrons in elastic collisions 
with neutral hydrogen atoms per unit mass is 



ricrn c T a T c . o. 
^ — ni{<T ea v ) 



(36) 



where m e and mu are the mass of electron and the mass 
of hydrogen atom, respectively, and the elastic scattering 
cross section is (Narita 1973) 



oo 



(cTeaV ) = / a ea v f (v) dv 



= 47ra; 



m 1/2 fkTA 3/2 r 



24 



(37) 



Here clq is the Bohr radius. 



Rate equations (23) and (|24j) imply that the number 
density of free electrons n e and atomic level populations 
rii change due to bound-free transitions, that is, due to 
ionizations and recombinations. The total ionization rate 



is 



Pik — n c Cik + i?ik , (38) 
where the rate of collisional ionizations is given by 



1/2 



T l/2 



exp 



Xi 



r, (T c 



(39) 



Xi = Xh/* 2 is energy of ionization from the i-ih level, 
Ti (T) is a slowly varying function of T evaluated with 
approximation formulae by Mihalas (1967). 
The rate of photoionizations is 



R t k = 4?r 



hp 



J v dv . 



(40) 



where abf {v, i) is an absorption cross-section at frequency 
v in bound-free transition from the z-th state and is a 
threshold frequency for ionization from the «-th state. 
The total recombination rate is 

Pki = n e C u + R ki , (41) 
where the collisional recombination rate is given by 



(33) Cm 



:Cik ■ 



The radiative recombination rate is 



Rki 



where 



-R 



Rl 



47T 



ik ! 



au\y, «) 



(42) 



(43) 



■ exp 



hv \ {2hv 3 , , , 



Substituting Eqs. (|42|) and (J43[) into Eq. 
that the total recombination rate is 



Pi 



|) we obtain 
(45) 



The system of ordinary differential equations ( p0| ) - 
( p4| ) is stiff because it is characterized by very different 
time constants due to the rate equations (23) and (pi|). 
In order to obtain the stable and enough correct solution 
of Eqs. ( p0| ) - (24) we used the Livermore solver for or- 
dinary differential equations based on the GEAR package 
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4. The radiation transfer equation 

For the shock wave models considered in the present study 
the extinction coefficient is highest at the Lyman edge 
frequency vqi. In the preshock region the corresponding 
extinction coefficient is xi^oi) ~ 3 ■ 10~ 4 cm" 1 . Behind 
the discontinuous jump the extinction coefficient initially 
decreases because of ionization of hydrogen atoms and 
then increases within the recombination zone. For all the 
models considered x(foi) & 5 ■ 10~ 3 cm -1 . Thus, the 
time needed for photons to travel a mean free path is 

tph. = [cxC^oi)] -1 ~ 7 ' 10~ 9 s ; where c is the velocity of 
light. Because a characteristic structural length for Lyman 
photons is only a few mean free paths, a photon flight time 
is much shorter than the characteristic time of producing 
the changes due to hydrodynamic motions and, therefore, 
the time-derivative term in the transfer equation can be 
omitted because the radiation field is quasi static to a 
good accuracy. 

In plane-parallel geometry the quasi static radiative 
transfer equation is written as 



S v . 



(46) 



where I v is the specific intensity of radiation depending on 
the frequency v, the directional cosine fi and the mono- 
chromatic optical depth dr v — \„dX. 

The present study is confined by treatment of the con- 
tinuum radiation transfer with total extinction coefficient 
given by 



(47) 



where Kbf(y) and K$(v) are bound-free and free- free ab- 
sorption coefficients, respectively. They can be found, for 
instance, in Mihalas (1978). or = 6.65 • 10 -25 cm 2 is 
the Thomson scattering cross section and ctr(^) is the 
Rayleigh scattering cross section evaluated due to approx- 
imate formulae by Kurucz (1970). Because the scattering 
is assumed to be isotropic, the source function S v can be 
written as 



Si/ 



r\v_ 

Xv 



(48) 



The source function S v is determined from integration of 
the system of ordinary differential equations ( p(i| ) - (£24J) 
and is evaluated at each cell center Xj_ 1 / 2 as a function 
of frequency v. 

Solution of the radiation transfer equation when ap- 
plied to the radiative shock waves is accompanied by some 
difficulties. First, because the shock wave propagates in 
the nearly neutral hydrogen gas, the slab has the apprecia- 
bly large optical depth in the Lyman continuum (r ~ 10 2 ), 
whereas in the Balmer and higher order continua the to- 
tal optical depth is very small. For example, the ratio of 
the total optical depth at the Lyman continuum edge to 



in order to obtain the stable solution for the entire fre- 
quency range we have to treat the transfer equation as 
a two-point boundary value problem. Second, the optical 
depth increments 

Atvj = i (w„ j _i/ 2 Ara j _i/2 + uj vj+1 / 2 Am j+1 / 2 ) , (49) 



are extremely small for hydrogen continua of order i > 3. 
Here w„ = Xvj P and Am^!^ = Pj-\i 2 AXj_x /% is the 
column mass contained in the cell. The straightforward 
application of the Feautrier method fails because of the 
limited machine accuracy. To alleviate this obstacle we 
employ an improved Feautrier solution proposed by Ry- 
bicki & Hummer (1991) and providing with much bet- 
ter numerical conditioning of the recurrence elimination 
scheme. 

Thus, the transfer equation ( ^6| ) is transformed (see, 
for example, Mihalas 1978) into the second-order differen- 
tial equation 



d 2 u 



Si/ ■ 



where u^ v is a mean-intensity-like variable defined by 



(50) 



' / 1 V 



= l[I(ji,u)+I(-n,u)] , 



(51) 



and p, changes in the range < p < 1. 

The transfer equation ( |50| ) is subject to boundary con- 
ditions at both surfaces of the slab. Assuming that radia- 
tion is produced only by the shock wave and that there is 
no incoming radiation from external sources, we have 



/„(//) = forX = Xx, 
I u l-H)=Q forX=X N , 



(52) 



The boundary conditions (|5j) are used in the second- 
order accuracy Taylor's expansion of the mean-intensity- 
like variable u v at both boundaries of the slab and are 
sufficient to complete the system of the finite-difference 
transfer equations. 



The transfer equation (50) is solved each cycle of 
global iterations for mean-intensity-like variable u^ v de- 
fined at the cell centers. The total number of cells is 
500 < N < 1200 depending on the shock wave model. The 
frequency range Vb < v < va is divided into Nq intervals, 
where the upper boundary of the range is va = 10 16 Hz. 
Boundaries of the intervals correspond to the threshold 
ionization frequencies and the lower boundary of the fre- 
quency range vb is the ionization threshold frequency of 
the A^c-th bound state. Within each interval the integral 
with respect to frequency v is replaced by the Gaussian 
quadrature sum, so that the integral over the whole fre- 
quency range [i/g, va\ is obtained by summation of interval 
integrals. The angular range < p < 1 is also replaced 
by a set of angular points {pi} at the Gaussian quadra- 
ture nodes. In the present study the number of quadrature 
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angular points were Nq = N m — 4. Integrating u v with 
respect to fi and v we evaluated in each cell center the 
mean intensity J„, the radiation energy density _E ra d, the 
radiation pressure P ra d and the divergence of radiative 
flux V-F r . The total radiative flux was calculated in cell 
boundaries from 



An J J v^vfi dfi dv . 
o o 



where 



(53) 



(54) 



is an antisymmetric flux-like variable. It should be noted 
that according to our conventions the radiative flux is 
negative if the radiation propagates in negative direction. 
Thus, the radiative flux emerging ahead the discontinuous 
jump is always negative. 

5. The global iteration procedure 

For integration of the system of ordinary differentail equa- 
tions (^) - ( pi|) we have to know the mean intensity J„ 
and the divergence of radiation flux V- F r at each cell cen- 
ter of the shock wave model. On the other hand, the radi- 
ation transfer equation (^0|) can be solved only when the 
extinction and emission coefficients are given in each cell 
center. In order to take into account a coupling between 
gas material and radiation field we employ the iteration 
procedure. If the initial approximation is enough close to 
the final solution, we may hope that each iteration will 
give a better approximation for the final solution than the 
previous one. In Fig. |l| is shown the flow chart where the 
main steps of the global iteration procedure are depicted. 

Within the whole shock wave both the radiation field 
and atomic level populations are in a strong departure 
from LTE. This feature is the principal difficulty accom- 
panying the shock wave model calculations because it is 
responsible for the narrow convergence area of global it- 
erations. In particular, the use of the initial LTE approxi- 
mation allows the converged solution to be obtained only 
for the weak shock waves with upstream Mach numbers 
Mi 3. For larger upstream velocities the initial oscilla- 
tion amplitude of the solution vector Y becomes so large 
that some quantities fall beyond their physical meaning. 
In order to alleviate this difficulty and to be able to con- 
sider the structure of stronger shock waves we computed 
a grid of the shock wave models with gradually increasing 
upstream Mach number, the LTE initial approximation 
being used only for the first model with the Mach number 
Mi = 2.3. Thus, each model of the grid with exception 
of the first one was computed with initial approximation 
obtained from the previous converged model having some- 
what smaller upstream velocity. The upstream velocity in- 



input of 
model parameters: 

Ui, T u Pl 



initialization of zonal quantities: 
LTE approximation for the 1-st model 
read the previous model otherwise 



Radiation Transfer 



compute of 

Xv, K. v , &i>, Si,, ■ 



solution of the 
transfer equation 
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ODE integration 
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Rankine-Hugoniot 
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ODE integration 







global iterations 
converged 





Yes 


output of the 
shock wave model 



Fig. 1. The flow chart of the global iteration procedure 

For test purposes some models were computed with differ- 
ent initial conditions taken from shock wave models with 
both larger and smaller upstream velocities. In all such 
cases the shock wave structure was found to converge to 
s single solution. 

In order to describe the convergence of global iterations 
we introduce for each component of the solution vector Y 
at the ^-th iteration the quantity 



e — max 

l<j<N 



1 



Vi 



specifying the maximum absolute relative change of the 
variable y across the whole shock wave model. A typical 
behaviour of global iterations is shown in Fig. || displaying 
in the upper panel the maximum absolute relative change 
of the electron temperature and the divergence of radia- 
tive flux within the whole shock wave model. On the lower 
panel of Fig. || is shown the iterative change of the total 
radiative flux F T emerging ahead the shock wave. Thus, 
if the initial approximation is good enough, we obtained 
in all cell centers of the model the exponentially decaying 
oscillations of the solution vector. The global iterations 



8 



Yu.A. Fadeyev & D. Gillet: The Structure of Radiative Shock Waves. I 



1 


] 

: i, 
: !„i 













- iliji 
: >;u 


r t 

;i i 
'i - ■ 








CO 












en 
_o 

-1 








r ■ * 

, . ■/ i 








\ v v J 






-2 


- 








i - 







50 






100 



E - 

U X 



I 




50 
iterations 



100 



Fig. 2. Upper panel: convergence plots for the electron 
temperature (solid line) and the divergence of radiative 
flux (dot-dashed line) given by the maximum absolute 
relative change e. Lower panel: convergence plot for the 
total radiative flux emerging from the outer boundary of 
the preshock region 

decrease. The final relative changes of the solution vec- 
tor depend on the both spatial resolution of the discrete 
model and accuracy of integration of differential equations 
( pfj| ) - (p4|). It should be noted also that the convergence 
of global iterations appreciably varies for different regions 
of the shock wave. In particular, the convergence is best in 
the radiative precursor and thermalization zone, whereas 
in the recombination zone becomes slower. Thus, the con- 
vergence plots shown in Fig. || display mostly variations of 
the solution vector in the vicinity of the outer boundary 
of the postshock region. 

6. Results of calculations 

In this paper we discuss the results of calculations done 
for the shock waves with upstream velocities 15 km s _1 < 
U\ < 60 km s _1 propagating through the unpertubed 




2x10" 



Fig. 3. The divergence of radiative flux against the spatial 
coordinate in the shock wave model with U\ = 40 km s _1 



pi = 10~ 10 gm cm" 3 (n H = 6.02 ■ 10 13 cm" 3 ). In to- 
tal we computed 46 models with upstream velocity incre- 
ment of AUi = 1 km s _1 . The outer boundary of the 
preshock region, where the gas is assumed to be unper- 
turbed, is set at X w —9.2 • 10 5 cm. Calculations were 
done for the two-level hydrogen atom, the first atomic 
state being treated in non-LTE. The radiation transfer 
equation was solved for the both Lyman and Balmer con- 
tinua (Nq — 2). Thus, the frequency point nearest to the 
Lyman edge frequency vqi = 3.288 • 10 15 Hz was set at 
i^q X = 3.754 ■ 10 15 Hz. More extensive calculations for the 
larger number of bound atomic states L and hydrogen 
continua Nq as well as for various temperatures T\ and 
densities p\ of the unperturbed gas will be given in the 
forthcoming paper. 

The most fascinating feature of radiative shock waves 
is that they demonstrate the strong interaction between 
gas material flows and the radiation field which they pro- 
duce. This interplay is best seen from the plots of the 
divergence of radiative flux as a function of the spatial 
coordinate X. One of such plots is shown in Fig. ||. By 
definition, the divergence of radiative flux is negative, if 
the fluid absorbs more energy than emits and, therefore, 
is heated. And conversely, when V-F r > 0, the gas radia- 
tively cools since it radiates more energy than it absorbs. 

As is seen in Fig. ||, the divergence of radiative flux 
is always negative in the preshock zone, the departure 
of V-F r from zero gradually increasing while the gas ap- 
proaches the discontinuous jump. Heating of the precursor 
gas material is due to absorption of the Lyman continuum 
radiation, hence the region, where the divergence of radia- 
tive flux V-F r perceptibly deviates from zero, extends over 
a few units of the optical depth at frequency v' Q1 . 

The properties of the radiation field does not change 
appreciably across the discontinuous jump. The spatial 
resolution of our shock wave models near the discontinu- 
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0.5 cm. The change of F T and -V-F r over this interval in- 
creases with increasing upstream velocity U± but does not 
exceed 0.3% for Ui = 60 km s _1 . Thus, just behind the 
discontinuous jump the divergence of radiative flux goes 
on to gradually decrease and reaches the minimum in the 
vicinity of the maximum of the electron temperature T e . 
The rapid growth of the divergence of radiative flux be- 
hind its minimum implies that the gas material flows into 
the radiatively cooling zone. For the shock wave model 
with U% = 40 km s _1 the distance between maximum and 
minimum of V-F r is ks 2.2 ■ 10 3 cm and corresponds to 
the optical depth between these layers of t{v' 01 ) ~ 3.5. It 
should be noted that both the minimum and maximum of 
the divergence of radiative flux as well as other properties 
of the shock wave very strongly (for some variables nearly 
exponentially) depend on the upstream velocity U\. 

6.1. The radiative precursor 

At the outer boundary of the preshock region the first 
level of the hydrogen atom is approximately in LTE since 
the Lyman continuum radiation emerging from the post- 
shock region is negligible. Both the Saha-Boltzmann re- 
lation and equations of statistical equilibrium give nearly 
the same (within a few per cent) the hydrogen ioniza- 
tion degree of xn ~ 10~ 8 . The main sources of opacity 
are bound-free transitions in the Lyman continuum (i.e. 
at frequencies v > uq\ = 3.288 • 10 15 Hz) and at lower 
frequencies the Rayleigh scattering by hydrogen atoms in 
the ground state. The free-free opacity and the Thomson 
scattering are negligible. With approaching to the discon- 
tinuous jump the hydrogen ionization degree increases, so 
that the both free-free opacity and Thomson scattering 
increase but nevertheless even for U\ = 60 km s" 1 they 
remain negligible within the entire preshock region. 

The radiative precursor is revealed as the part of the 
preshock region where the hydrogen gas is heated and is 
ionized by the Lyman continuum radiation emerging from 
the postshock region. The temperature of heavy particles 
(neutral hydrogen atoms and hydrogen ions) T a remains 
constant (T a = T\ = 3000K) in the radiative precursor. 
In Fig. ^ are shown the plots of the electron temperature 
T e and the hydrogen ionization degree xn as a function 
of the distance from the discontinuous jump for the shock 
wave models with upstream velocities U% — 40, 45, 50 
and 55 km s _1 . As is seen, the perceptible heating and 
ionization occur at distances from the discontinuous jump 
smaller than X prcc w 3 • 10 4 cm. The distance X prcc cor- 
responds to the optical depth of t(u' q1 ) ~ 3. Thus, for 
Ui < 60 km s _1 the geometrical thickness of the radiative 
precursor approximately does not depend on the upstream 
velocity. 

Within the radiative precursor the hydrogen gas and 
radiation field are in strong departure from LTE. In partic- 
ular, ionization of hydrogen atoms is mainly due to radia- 




2 4 

og X (cm) 



Fig. 4. The hydrogen ionization degree cch (upper panel) 
and the electron temperature T c (lower panel) in the 
preshock region of the shock wave models with upstream 
velocities U\ = 40, 45, 50 and 55 km s _1 



ciably exceed collisional ionizations. Fig. |5| shows the rates 
of radiative transitions i?ik and i?ki = ( n i/ n o)^ik as weu 
as the collisional rates n e Cik and n e Cki = n e (n 1 /n*)Cik 
as a function of X for the shock wave model with U\ = 
40 km s _1 . The precursor transition rates are very sen- 
sitive to the upstream velocity XI \, For example, for up- 
stream velocities 20 km s _1 < U\ < 60 km s _1 the rates 
of ionizations from the ground state increase in the ranges 
-7.1 < logi?ik < 2.6 and -25.8 < log (n e Ci k ) < -10.6, 
respectively. However, notwithstanding such a strong de- 
pendence on Ui, for all shock wave models the spatial 
dependencies of transition rates were found to be qualita- 
tively similar to those displayed in Fig. 

The growth of the hydrogen ionization degree with ap- 
proaching to the discontinuous jump is accompanied by 
the increase of the gas pressure gradient. The correspond- 
ing decrease of the gas material velocity becomes, how- 
ever, perceptible only for U > 50 km s . For example, at 



10 



Yu.A. Fadeyev & D. Gillet: The Structure of Radiative Shock Waves. I 



o 



1 ' 

" R ,k 


" R k1 




" n a C k1 


\ \. 




\ 

\ \ 




\ \ - 




\ \ . 




\ _ _ _ _ 


- n e C 1k 




I , I , l , I - 



J i I i L 

2 4 

log X (cm) 



Fig. 5. The rates of ionizations (solid lines) and recombi- 
nations (dashed lines) against the distance from the dis- 
continuous jump in the pre-shock region of the shock wave 
model with U\ = 40 km s _1 

Table 1. The properties of the preshock inner boundary 



log (transition rate s j 
Ui logXH To i?ik -Rki n c Cik n Cki 



15 


-8.04 


3000 


-7.36 


-8.79 


-25.85 


-17.86 


20 


-8.04 


3000 


-7.11 


-8.79 


-25.85 


-17.86 


25 


-8.04 


3000 


-7.01 


-8.79 


-25.85 


-17.86 


30 


-7.21 


3000 


-4.86 


-7.99 


-25.05 


-16.26 


35 


-4.19 


3002 


-1.20 


-4.93 


-21.98 


-10.16 


40 


-2.92 


3018 


0.23 


-3.65 


-20.58 


-7.61 


45 


-2.33 


3062 


0.90 


-3.03 


-19.67 


-6.45 


50 


-1.91 


3158 


1.40 


-2.53 


-18.56 


-5.62 


55 


-1.47 


3416 


1.90 


-1.91 


-16.46 


-4.77 


60 


-0.84 


4590 


2.63 


-0.75 


-10.61 


-3.62 



crease of U is w 0.15% for U\ — 55 km s 1 and is « 0.5% 
for Ui = 60 km s _1 . 

Of great interest are the physical conditions at the in- 
ner boundary of the preshock region. For the models con- 
sidered in the present study the spatial coordinate of this 
boundary is Xj_i/2 = —0.5 cm. For upstream velocities 
U\ SS 28 km s _1 both the hydrogen ionization degree and 
the electron temperature nearly do not change within the 
preshock region. These quantities show the perceptible de- 
pendence on upstream velocity only for U\ > 28 km s _1 . 
In Tabled we give the hydrogen ionization degree ih, the 
electron temperature T e as well as the transition rates at 
the inner boundary of the preshock region. For the sake of 
convenience the upstream velocity U\ is given in km s _1 . 

The electron temperature just ahead the discontinuous 
jump increases with increasing upstream velocity nearly 
exponentially. For the shock wave models of the present 



due to the following fitting formula 

logT c = 3.303 + 2.115 • 10~ 2 [/i - 7.590 • lO" 4 ^ 2 + 

+ 8.326 • lO" 6 ?/ 3 , (55) 

where the upstream velocity U\ is expressed, for the sake 
of convenience, in km s _1 . 

As is seen in Table [j], the dominating process just 
ahead the discontinuous jump is photoionization and, the- 
refore, the shortest relaxation time in the radiative pre- 
cursor is that of photoionizations from the ground state 
£phi = 1/-Rik- The photoionization relaxation time grad- 
ually decreases with increasing upstream velocity from 
iphi « 2 • 10 7 s for [/i = 15 km s" 1 to t phi w 2 • 10~ 3 s 
for Ui = 60 km s _1 . Comparing these relaxation times 
with the time needed for gas to flow through the precur- 
sor ihyd = X picc /Ui we find that for models with upstream 
velocities U\ < 50 km s _1 the ratio of the photoionization 
relaxation time to the hydrodynamic time is tphiAhyd ^ 1 
and only for U\ = 60 km s _1 the photoionization relax- 
ation time becomes nearly comparable with hydrodynamic 
time: iphiAhyd ~ 0.4. Because for establishment of the 
statistical equilibrium this ratio should be tphiAhyd 1, 
the ground state populations of the hydrogen atom are 
in strong departure from the statistical equilibrium and 
the hydrogen ionization degree cannot be described in as- 
sumption of statistical equilibrium. For example, at the 
inner boundary of the preshock region of the shock wave 
model with U\ = 50 km s _1 the hydrogen ionization de- 
gree is in = 1.22 • 10 -2 , whereas solution of the equations 
of statistical equilibrium gives xh = 0.99. 

6.2. The thermalization and recombination zones 

Behind the discontinuous jump the translational kinetic 
energy of heavy particles is redistributed among various 
degrees of freedom characterized by different relaxation 
times. The fastest relaxation process is the translational 
kinetic energy exchange in clastic collisions of electrons 
with neutral atoms and ions. Another relaxation process is 
excitation of bound atomic states and ionization of hydro- 
gen atoms. Both excitation and ionization need, however, 
much more collisions than translational kinetic energy ex- 
change (Stupochenko et al. 1967), so that just behind 
the discontinuous jump the electron temperature grad- 
ually increases, whereas the hydrogen ionization degree 
remains nearly constant. Note that although the bound- 
bound transitions were not considered in our model, the 
excitation of atomic states is taken into account as a result 
of bound-free transitions. 

In Fig. ^ are shown the electron temperature T c and 
the temperature of heavy particles T a as a function of dis- 
tance from the discontinuous jump for shock wave models 
with upstream velocities from 15 km s _1 to 60 km s _1 . 
As is seen, the characteristic time of the electron temper- 
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Fig. 6. The electron temperature T (solid lines) and tem- 
perature of heavy particles T a (dashed lines) as a function 
of distance from the discontinuous jump in the postshock 
region 

velocity. Furthermore, for U± < 30 km s _1 there is a tem- 
perature plateau within of which both temperatures ap- 
proximately do not change until they begin to decrease. 
The existence of the temperature plateau is due to the fact 
that for upstream velocities U\ < 30 km s _1 the electrons 
acquire the energy from heavy particles due to elastic col- 
lisions with neutral hydrogen atoms. The rate of energy 
gain in such collisions is very small and gradually decreases 
when the electron temperature T c approaches the temper- 
ature of heavy particles T a . The temperature plateau ap- 
pears when the electron energy gain in elastic collisions 
with neutral hydrogen atoms becomes almost negligible. 
At upstream velocities U\ > 30 km s^ 1 the electron tem- 
perature plateau is ended by the slight bump, the bump 
being wider and higher with increasing U±. 

In the second and third columns of Table || are given 
both the maximum value of the electron temperature in 
the postshock region as well as the distance of this point 
from the discontiuous jump expressed in cm. In the first 



Ui 


Tc 


logX 


xn 


logX 


log±V-F r 




15 


7521 




0.014 


6.29 


11.097 


-5.22 


20 


11486 




0.071 


5.72 


12.671 


-3.43 


25 


16281 


4.95 


0.134 


5.16 


13.201 


-2.83 


30 


20986 


4.20 


0.206 


4.40 


13.473 


-2.80 


35 


24712 


3.42 


0.294 


3.89 


13.749 


-2.19 


40 


28274 


2.92 


0.406 


3.61 


13.992 


-1.81 


45 


31681 


2.58 


0.528 


3.43 


14.198 


-1.66 


50 


35050 


2.29 


0.662 


3.32 


14.381 


-1.49 


55 


38575 


1.96 


0.805 


3.22 


14.569 


-1.23 


60 


43371 


1.40 


0.964 


2.89 


14.848 


-0.74 



for the sake of convenience, in km s _1 . The electron tem- 
perature peak was not found for the shock wave mod- 
els with U\ < 20 km s _1 and in these cases we give 
only the electron temperature of the plateau. The elec- 
tron temperature maximum can be approximately con- 
sidered as a point where temperatures of heavy particles 
and free electrons equilize. Thus, the width of the relax- 
ation zone where both temperatures equalize decreases by 
a factor of » 3600 for the upstream velocity increasing 
from 25 km s _1 to 60 km s _1 . It is of interest to note that 
for upstream velocities U\ < 50 km s _1 the temperature of 
heavy particles T a remains nearly constant until the elec- 
tron temperature begins to increase just before its drop. 
This is due to the fact that at small upstream velocities 
the fractional abundance of free electrons is so small that 
they cannot perceptibly affect the gas of heavy particles. 

The translational energy exchange between heavy par- 
ticles and electrons is due to elastic collisions of electrons 
with both neutral hydrogen atoms and hydrogen ions. 
The cross section of elastic collisions with hydrogen ions 
is much larger than that of electrons with neutral atoms 
and, therefore, the translational energy gain by electrons 
from heavy particles is strongly dependent on the hydro- 
gen ionization degree. The rates of energy gain by elec- 
trons in both these processes are shown in Fig. [?]. The 
abrupt decrease of Q C j and Q ea occurs when both temper- 
atures equalize. 

For shock wave models considered in the present study 
the rate of energy gain by electrons in elastic collisions 
with ions exceeds that in elastic collisions with neu- 
tral atroms only when the hydrogen ionization degree is 
ih > 10~ 2 . As is seen in Fig. [?], for upstream velocities 
Ui < 25 km s _1 , the electron temperature equalizes with 
temperature of heavy particles only due to elastic colli- 
sions of electrons with neutral atoms. Comparing with 
lower panel of Fig. ^| one sees that when the electron tem- 
perature reaches the plateau, the rate of energy gain Q ea 
decreases by nearly two orders of magnitude. 

At upstream velocities 30 km s _1 < U\ < 50 km s _1 , 
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Fig. 7. The rates of energy gain by electrons in elastic 
collisions with hydrogen ions (solid lines) and neutral hy- 
drogen atoms (dashed lines) 



dominate only just before the temperature drop. It is 
this increase of Q el that is responsible for the electron 
temperature peak near the end of the electron tempera- 
ture plateau. The translational energy gain by electrons in 
elastic collisions with hydrogen ions becomes completely 
dominating for U\ > 50 km s _1 . The gradual decrease of 
Q cl just behind the discontinuous jump in the shock wave 
models with Ui > 40 km s _1 (see upper panel of Fig. 
is due to the temperature dependence of the equipartition 
time given in Eq. (34). 

In Fig. 8 are shown the plots of the hydrogen ionization 
degree in the postshock region of the shock wave mod- 
els with upstream velocities U\ — 40, 50 and 60 km s _1 . 
Comparing with Fig. ^ one sees that the maximum of the 
hydrogen ionization occurs at much larger distances from 
the discontinuous jump than the maximum of the electron 
temperature. Very approximately the distances of both 
these maxima can be expressed as power functions of the 
upstream velocity: 



"0.5 




og X (cm) 

Fig. 8. The hydrogen ionization degree as a function of 
the distance from the discontinuous jump in the shock 
wave models with upstream velocities U\ = 40, 50 and 
60 km s _1 



X(maxx H ) = 6.07- lO 7 ^ 



71-7—5.60 



(57) 



In these approximate expressions, for the sake of conve- 
nience, the distances AT(maxT c ) and X(maxxn) are ex- 
pressed in km whereas the upstream velocity is expressed 
in km s . The maximum values of xh as well as the corre- 
sponding distances from the discontinuous jump expressed 
in cm are given in the fourth and fifth colums of Table ^, 
respectively. According to the schematical division of the 
shock wave noted above the maximum of can be con- 
sidered as the boundary between the thermalization and 
recombination zones. 

The large distance between maxima of the electron 
temperature and the hydrogen ionization degree implies 
that the degree of freedom associated with ionization of 
hydrogen atoms is frozen in comparison with excitation of 
translational motions. As a result, the gas flows through 
the maximum of T e at smaller heat capacity C p than it 
would be in equilibrium. Thus, when the perceptible frac- 
tion of hydrogen atoms is ionized, the gas density begins 
to increase. At larger distances the gas density goes on to 
increase due to the radiative cooling (see Fig. ^|). 

The total ionization rate of hydrogen atoms is a sum 
of photoionizations and collisional ionizations, so that it 
depends on the number density of free electrons. Because 
no ionization occurs across the discontinuous jump, the to- 
tal ionization rate behind the discontinuous jump strongly 
depends on the ionization in the radiative precursor. To 
demonstrate this dependence, in Fig. [l^ are shown the 
rates of ionizations and recombinations, both collisional 
and radiative, in the postshock regions of the shock waves 
with upstream velocities 30 and 50 km s . 

For the shock wave model with U\ = 30 km s^ 1 the 
hydrogen ionization degree at the discontinuous jump is 
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Fig. 9. The postshock gas density p against the distance 
from the discontinuous jump. The numbers at the curves 
indicate the upstream velocity Ui in km s _1 



the role of collisional ionizations is negligible and the num- 
ber density of free electrons increases initially only due to 
photoionizations. For U\ = 50 km s _1 the hydrogen ion- 
ization degree at the discontinuous jump is ih ~ 1CP 2 . 
Free electrons created in the precursor play a role of the 
seeds producing yet more electrons and leading to the elec- 
tron avalanche. 

Compared to the preshock region where the contri- 
bution of opacity sources nearly does not depend on the 
upstream velocity, the gradual growth of the hydrogen 
ionization degree behind the discontinuous jump leads to 
the appreciable dependence of postshock opacities on the 
upstream velocity U\. For XJ\ < 30 km s _1 the post- 
shock opacity in the Balmer continuum is mainly due 
to Rayleigh scattering by neutral hydrogen atoms. Only 
near the Balmer edge, the opacity due to bound-free tran- 
sitions from the second level becomes most important. 
With increasing upstream velocity the role of the bound- 
free transitions from the second level increases and for 
U\ > 50 km s _1 this opacity mechanism becomes domi- 
nating within the whole range of the Balmer continuum. 
The role of the Thomson scattering is somewhat percepti- 
ble only near the Balmer edge, whereas the free-free opac- 
ity can be neglected. 

6. 3. The radiation field 

Although we did not include into the system of ordinary 
differential equations (|2C| ) - (^4|) the terms with radia- 
tion energy density i? ra( j and radiation pressure P ra d, these 
quantities were evaluated as 

oo 1 

£rad = — / / u^dpdv (58) 




2 




og X (cm) 

Fig. 10. The rates of ionizations (solid lines) and recom- 
binations (dashed lines) against the distance from the dis- 
continuous jump X in the postshock region. Upper panel: 
U\ = 50 km s _1 ; lower panel: U\ = 30 km s _1 

oo 1 

frad = -^7 J J u^^dp.dv (59) 



together with total radiative flux F r and divergence of the 
radiative flux V-F r each time when the radiation trans- 
fer equation was solved. For the fixed spatial coordinate 
X the ratio of the radiation energy density to the total 
translational kinetic energy gradually increases with in- 
creasing upstream velocity and is in the range 5 • 10~ 4 Si 
E I&d /{E a + E c ) <, 0.02 for 20 km s" 1 < Ui < 60 km s _1 . 
The ratio E Ya d/(E a + E C ) is highest at the inner boundary 
of the preshock region because just behind the discontin- 
uous jump the total translational kinetic energy increases 
by more than an order of magnitude, whereas the total 
change of E la ^ within the shock wave does not exceed 
30%. Because the most of the energy flux is contained in 
the radiative flux F T , we find that UE lad /F r < 2.3 • 10~ 4 
for Ui — 60 km s _1 . Thus, our assumption that the ra- 
diation energy and radiation pressure can be neglected is 
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log X (cm) 

Fig. 11. The total radiation energy density as a function 
of distance from the discontinuous jump in the shock wave 
model with U\ = 60 km s _1 

In Fig. [ll] are shown the dependencies of -E ra d on the 
distance from the discontinuous jump for the both pre- 
shock and postshock regions of the model with U\ = 
60 km s _1 . The increase of -E ra d in the preshock region 
with approaching to the discontinuous jump is due to the 
radiative heating by the layers absorbing the Lyman con- 
tinuum radiation. Behind the discontinuous jump the ra- 
diation energy density reaches the maximum in the hy- 
drogen recombination zone. 

The most of the shock wave radiation, by definition, 
is produced in the layers where the divergence of radia- 
tive flux V-F r reaches the maximum. These layers locate 
nearly at the same distance from the discontinuous jump 
as the maximum of the hydrogen ionization degree. The 
maximum of iV-F r exponentially increases with increas- 
ing upstream velocity and for 20 km s _1 < U\ < 60 km s" 1 
is given by 

(-V-F r ) = 1.25-10 12 exp(0.105£/i) erg gm~ V\(60) 

\ r / max 

where the upstream velocity is expressed in km s . Max- 
imum values of -V-F r are given in the sixth column of 
Table|| 

As was noted above the shock wave models have a 
small optical depth in the Balmer continuum (T max ~ 
10~ 4 ) and at the same time they are opaque for the Ly- 
man continuum photons. Therefore, the radiation emerg- 
ing from the both surfaces of the slab is mostly within 
the Balmer continuum whereas the Lyman continuum ra- 
diation is transported only within the narrow zone sur- 
rounding the discontinuous jump. In the lower panel of 
Fig. [[2] for the shock wave model with upstream velocity 
Ui = 40 km s -1 is shown the Lyman continuum radiative 
flux Fl y as a function of the spatial coordinate X. This de- 




X (cm) 

Fig. 12. The radiative flux within the Lyman continuum 
i<Ly (lower panel) and the total radiative flux F r (upper 
panel) against the spatial coordinate X for the shock wave 
model with U\ = 40 km s _1 



of radiative flux shown in Fig. The negative flux corre- 
sponds to the radiation propagating in the direction of the 
preshock region. The optical depth between minimum and 
maximum of Fl y is IS.t(v' 01 ) « 4.4 for U\ — 40 km s _1 and 
slowly decreases from At(v' 01 ) « 24 for U\ = 30 km s _1 to 
At(i/' 01 ) w 1 for U\ — 60 km s _1 . This decrease is mostly 
due to the growth of the hydrogen ionization degree in 
the postshock region. The ratios of the maximum radia- 
tion flux transported within the Lyman continuum to the 
total radiation flux are given in the last column of Table |[ 
In the upper panel of Fig. [l^ for the same model is 
shown the total radiative flux F T . This plot displays only 
the very vicinity of the discontinuous jump. In the pre- 
shock region, beyond the radiative precursor, the total ra- 
diative flux F r remains constant since the divergence of 
radiative flux is V-F r = 0. In the postshock region the 
total radiative flux gradually increases and asymtotically 
tends to the same value as in the preshock but with the 
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the total radiative fluxes emerging from both surfaces of 
the slab are — F t i = F t n- 

7. Conclusion 

The primary goal of the work reported in this paper was to 
obtain the self-consistent solution of the equations of fluid 
dynamics, rate equations and radiation transfer equation 
for the structure of the steady shock wave. The procedure 
of global iterations described in the present paper in gen- 
eral resembles the compute of stellar atmosphere models. 
Indeed, like in stellar atmosphere calculations the shock 
wave model takes into account the coupling between the 
gas material and radiation field. The self-consistent model 
is obtained with iteration procedure comprising the solu- 
tion of the radiation transfer equation and integration of 
the mass, momentum and energy conservation equations 
written in the form of the ordinary differential equations. 
Each cycle of iterations gives, in general, improved char- 
acteristics of the gas and radiation field. 

At the same time, the problem of the shock wave struc- 
ture compared to that of stellar atmosphere models con- 
tains a number of serious complications. First, atomic level 
populations are not only in strong departures from LTE 
but are also in significant departures from statistical equi- 
librium. Second, unlike the stellar atmospheres, where the 
divergence of radiative flux is V- F r = (the condition of 
radiative equilibrium), in shock waves the part of the en- 
ergy of hydrodynamic flow is transformed into radiation 
and the radiative equilibrium is established only far away 
from the discontinuous jump. Furthemore, in stellar at- 
mosphere models the total radiative flux is given as one 
of the boundary conditions, whereas in the shock wave 
model the emerging flux is obtained from the solution of 
the problem. The small optical depth increments in hy- 
drogen continua of order I > 3 lead to the losses of the 
machine accuracy when the Feautrier technique is applied. 
The problem is so serious, that even the improved method 
by Rybicki & Hummer (1991) sometimes fails. Third, the 
rate equations are stiff and need the special treatment in 
their solution. In particular, the convergence of global it- 
erations depends on the tolerance parameter determining 
the maximum error permitted during the integration. 

The present paper is confined by consideration of the 
two-level atomic model, so that the radiation transfer is 
treated for the Lyman and Balmcr continua, only. This 
approximation seems to be insufficient for the shock wave 
problem because the occupation numbers of levels I > 2 
obviously deviate from LTE and the perceptible frac- 
tion of radiation is transported at frequencies lower than 
the Balmer edge frequency. Nevertheless, notwithstand- 
ing such a restriction, there is a qualitative agreement of 
our results with those obtained earlier by other authors. 
For example, according to calculations of Gillet & Lafon 
(1990) the electron temperature just ahead the discon- 



of Ui — 80 km s _1 . Although in the present study the 
highest upstream velocity was U\ = 60 km s _1 , a very 
approximate comparison can be done with fitting formula 
( |55|) which gives the same electron temperature for the 
upstream velocity U\ « 75 km s _1 . Thus, more detailed 
calculations are needed and in the forthcoming paper we 
are going to present the grid of the shock wave models 
computed for the larger number of hydrogen atomic levels 
and wider range of upstream velocities, the models of the 
present study being used as initial approximation for more 
correct shock wave models. 

More realistic models, however, should take into ac- 
count not only bound-free terms but also bound-bound 
terms in the rate equations and the radiation transfer 
problem should be solved for the both continuum and 
spectral line radiation. This is the perspective for the near 
future. It is certainly one of the most basic. Indeed, pre- 
ceding shock studies show that radiative processes, which 
determine the wake cooling, have a strong influence on 
the resulting shock structure. Because in the model of this 
paper we consider a pure H-plasma without H~ , and only 
include the bound-free photo- and collisional processes of 
H atoms, we expect that the absence of some predom- 
inant coolants such as neutral and singly ionized metal 
atoms might appreciably underestimate the radiative cool- 
ing rate of the gas. The importance of radiative heating 
and cooling rates in shocked circumstellar envelopes have 
been recently investigated (Woitke et al. 1996) but only a 
few transitions of the numerous metal lines were consid- 
ered. At present such a basic study seems to be beyond 
our immediate abilities and is out of the scope of our pa- 
per which was to investigate the possibility of obtaining a 
self-consistent solution of the structure of radiative shock 
waves in dense atmospheric gas. 
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